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c/3 \ Abstract 

We model a planar array of fluxlines as a discrete solid-on-solid model with quenched 
J5 ■ disorder. Simulations at finite temperatures are made possible by a new algorithm 

which circumvents the slow glassy dynamics encountered by traditional Metropolis 
n^ \ Monte Carlo algorithms. Numerical results on magnetic susceptibility variations 

Ch ' support analytic predictions. 
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(^ \ Over the years various low-temperature glass phases have been proposed for 

O ■ vortex matter in dirty type-II superconductors [1]. Understanding the thermo- 

Qv ■ dynamics of such interacting many-body systems in the presence of quenched 

Q\ . disorder remains a challenge. Similarities between the randomly-pinned vortex 



system and the more familiar mesoscopic electronic systems [2] are highlighted 

by a recent experiment on vortices threaded through a thin crystal of 2H-NbSe2 

by BoUe et al. [3]. Sample-dependent magnetic responses known as "magnetic 

2 ■ finger prints" have been observed for such a mesoscopic vortex system. 

O' 

^ ' The simplicity of the planar vortex system allowed detailed theoretical studies 

of this glassy system [4-9], with quantitative predictions including the uni- 
^ ■ versal sample-to-sample variation of magnetic responses. However, until the 

work of Bolle et al., there were hardly any experimental studies of this sys- 
tem, with difficulties stemming partly from the weak magnetic signals in such 
2d systems. Meanwhile, progress on numerical approaches were also seriously 
hampered by the slow glassy dynamics often encountered in simulations at 
finite temperatures [10]. To overcome the this numerical difficulty, we model 
a planar array of fluxlines as a discrete solid-on-solid model, viz., the dimer 
model with quenched disorder. A new polynomial algorithm circumventing the 
glassy dynamics enables us to simulate large systems of this discrete model at 
any flnite temperature. 
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Fig. 1. (a) Snapshot of a dinier covering (thick bonds) together with the associated 
"height" values ^(r) for a lattice of size L = 8 at temperature T^i = 1.0. (b) Dimer 
covering of the fixed reference, (c) Vortex line configuration (thick lines) obtained 
as the difference between (a) and (b). 

The Model: The dimer model consists of all complete dimer coverings {D} 
on a square lattice C as illustrated in Fig. 1(a). The partition function is 



^exp 

{D} 



J2 ^ij/'^d 
<ij>eD 



where the sum in the exponential is over all dimers of a given covering, and T^ 
is the dimer temperature. Quenched disorder is introduced via random bond 
energies e^j, chosen independently and uniformly in the interval 
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The dimer model is related to the planar vortex-line array via the well-known 
mapping to the solid-on-solid (SOS) model (see Figs. 1): Define integer heights 
{h} at the centers of squares of C which form the dual square lattice C^, then 
orient all bonds of C^ such that elementary squares of C^ that enclose sites of 
a chosen sublattice of C (solid dots in Fig. 1(a)) are circled counterclockwise; 
Assign —3 to the difference of neighboring heights along the oriented bonds 
if a dimer is crossed and -|-1 otherwise. This yields single-valued heights up 
to an overall constant. In terms of the height configuration h{r), the partition 
function (1) can be written alternatively as Z = J2{h(r)}^~^'^^'^\ where the 
SOS Hamiltonian takes the following form in the continuum limit, 
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{Vhf - f (r) -Vh + g cos {Gh) 



(2) 



Here, K is an effective stiffness caused by the inability of a tilted surface to 
take as much advantage of the low weight bonds as a flatter surface, and f (r) 
is a random local tilt bias. The periodicity of the cosine potential in (2) is 
given by G = 27r/4 since the smallest "step" of this height profile is four. 
In the present context of a randomly pinned vortex array, /i(r) describes the 
coarse-grained displacement field of the vortex array with respect to a uniform 
reference state at the same vortex line density; see Fig. 1(c) and Refs. [5,8]. 



The Algorithm: Computing the partition function of complete dimer cov- 
erings on a weighted p/anar lattice Q can be achieved in polynomial tirae [12]. 
Weights here refer to the Boltzmann factors Wij = exp{—eij/Ta) on the bonds. 
In this article, we shall use the shuffling algorithm by Propp and coworkers[13] 
to sample the configurational space. The algorithm relates the partition func- 
tion Zl on a lattice of linear size L to Zl^i as Zl{{w}) = C {{w}) Z l^i{{w'}) 
after a simple weight transformation {w} -^ {w'}. The prefactor C{{w}) is 
independent of dimer coverings. The partition function is obtained in a "de- 
flation" process in which the above recursive procedure is carried out down to 
L = with Zq = 1. This deflation process can also be reversed in an "infla- 
tion" process where a dimer covering at size L — 1 can be used to stochastically 
generate a dimer covering at size L according to Zl already obtained. Repeat- 
ing the inflation process thus generates uncorrelated "importance samplings" 
of the dimer conflgurations, or equivalently the equilibrium height conflgura- 
tions, without the need to run the slow relaxational dynamics. The ensuing 
numerical results are obtained by taking various measurements of the height 
conflgurations generated this way. A somewhat inconvenient feature of this ap- 
proach is that the algorithm requires open boundary condition on the dimer 
model; this in turn flxes the total number of vortex lines, e.g., to L/2 on a 
L X L lattice (see Fig. 1). 



Numerical results: Since the temperature of the vortex array, given by 
K^^ in (2), is generally different from the dimer temperature T^, we flrst need 
to calibrate the temperature scale. To do so, we exploit a statistical rota- 
tional symmetry [15,9] of the system (2) which guarantees that at large scales, 
the effective K is the same as that of the pure system. In particular, K can 
be obtained by measuring the disorder-averaged thermal fluctuation of the 

displacement held Wj. = (-< hj^{r) > < hxir) y"^) = {2irK)~^\n{L), where 

/iy(r) = h{r) — {h{r)) measures the thermal distortion superposed on the dis- 
torted "background" {h{r)) by the quenched disorder. We used -< ... y, (...), 
and overline to denote spatial, thermal, and disorder averages respectively. 
With the polynomial algorithm, we were able to perform thorough disorder 
averages for equilibrated systems of sizes up to 512 x 512. To reduce boundary 
effects, we focus on the central L/2 x L/2 piece of the system and compute its 
displacement fluctuations. Fig. 2(a) illustrates the dependence of W^{L) for 
various dimer temperatures T^. The linear dependence on ln(L) is apparent. 
Identifying the proportionality constant with {27rK)~^, we obtain the empir- 
ical relation between K~^ and T^ shown in Fig. 2(b). Note that our result 
recovers the exact relation K~^{Td — > oo) = IG/vr for the dimer model with- 
out disorder [14]. Since the glass transition of the system (2) is expected to 
occur at temperature K~^ = Air/G'^ = IG/tt, our system is glassy for the entire 
range of dimer temperatures. 

To probe this glassy state further, we now focus on its magnetic response, i.e.. 
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Fig. 2. (a) The disorder-averaged thermal displacement fluctuation W^(L) for var- 
ious dimer temperatures T^. Each data point was obtained by generating 10^ ther- 
mal samplings for each of 10'^ disorder realizations. The straight lines are fit to 
Wj^ = {27rK)~^ lii(-^) + const, (b) The extracted relation between K~^ and Td- The 
horizontal line indicates the asymptotic value of K~^ as Td 
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the uniform magnetic susceptibility x- Due to the constraint of fixed vortex 
density when using the dimer representation, it is not easy to compute the 
magnetic susceptibility directly by varying the vortex fugacity (i.e., external 
magnetic field). Instead, we use the fluctuation-dissipation relation. 



Xoclimq 



(/lq/l_q) - (/iq)(/i-q) 



limq^C2(q) 

q-*0 



(3) 



where h^ is the Fourier transform of h{r). x only probes the system at the 



largest scale as indicated by the q 
(vr/L, 0) for finite systems. 



limit, which is evaluated at (g^:, Qy) 



As shown in Figs. 3, the distribution P{x) at low temperature is peaked at zero 
and decays more strongly than any power laws for large x's. This indicates 
that the glass phase is typically rigid against the penetration of additional 
fiuxlines when the external magnetic field is increased by a small amount. 
The mean value x, which takes on a finite value oc K~^ due to the statistical 
rotational symmetry [9,15], is however controlled by the rare events contained 
in the tail of P{x)- This is similar to the behavior found for a single flux line 
in 1 + 1 dimensional random potential [16,17]. 

A striking feature of this vortex glass concerns the sample-to-sample variation 
of the susceptibility x- It was predicted [9] that the fractional variance. 



var[x]/x' 



.2 L- 



D\t\ for < -r < 1, 



(4) 



is universal, with D being a computable, size-independent constant of or- 
der unity. Using Eq. (3), we obtain var[x] oc \im.(^^Qq'^Ci{q) where C^i^q) = 

2 

CKq) — C'2(ci) • As shown in Fig. 4(a), q^C4(q) obeys a power-law scaling as 
q^C4(q) ~ (qL)^-'^^^'^^ for small q. This behavior is consistent with the analyti- 
cal prediction that the variance is size-independent because the q —^ limit is 
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Fig. 3. (a) q^C2iq) vs. q 



^ Different symbols denote different disorder realizations, 
all of which are for L = 128 at T^ = 0.04. The susceptibility x is evaluated at 
q = tt/L. (b) Probability distribution P{x) computed from 6 x 10^ disorder real- 
izations for L = 16 at Td = 0.04. Dotted line is a fit to a stretched exponential 
P{x) = ^exp(-Sx°) with a = 0.5(2). 
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Fig. 4. (a) g^L'(q) vs. {qLj^ along qy = Q axis. Data are obtained by averaging over 
10^ disorder realizations for L = 16, 32, 64, 128, 256 at T^ = 0.2. The solid line is the 
least square fit; (b) Fractional variance as a function of the reduced temperature, 
r = 1 — Kg/K^ for L = 32, 64, 128, and 256 averaged over 10^ disorder realizations. 
The solid line is a least-square linear fit for |r| < 0.15. 

evaluated at {qxiQy) = (tt/L, 0) for finite systems. The power-law scaling and 
its exponent are, however, not understood at present. Finally the temperature 
dependence of the fractional variance, as shown in Fig. 4(b), can be deduced 
from the relation between K~^ and T^ given in Fig. 2(b). It is well described 
by the linear form (4) close to the glass transition with D = 0.454(5). 



Conclusion. Finite temperature simulations were performed on a disor- 
dered dimer model to study the magnetic response of a planar vortex array. 
Universal susceptibility variations in the collective pinning regime were ob- 
served; critical behavior near the glass transition compared well with analytic 
predictions. It is hoped that the present study will stimulate further experi- 
mental investigations of the physics of mesoscopic vortex systems. 
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